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ABSTRACT. We present a high-resolution Gravity Recovery and Climate Experiment (GRACE) mascon 
solution for Gulf of Alaska (GOA) glaciers and compare this with in situ glaciological, climate and other 
remote-sensing observations. Our GRACE solution yields a GOA glacier mass balance of -65 ± 1 1 Gta -1 
for the period December 2003 to December 201 0, with summer balances driving the interannual 
variability. Between October/November 2003 and October 2009 we obtain a mass balance of 
-61±11Gta _1 from GRACE, which compares well with -65±12Gta _1 from ICESat based on 
hypsometric extrapolation of glacier elevation changes. We find that mean summer (June-August) air 
temperatures derived from both ground and lower-troposphere temperature records were good 
predictors of GRACE-derived summer mass balances, capturing 59% and 72% of the summer balance 
variability respectively. Large mass losses during 2009 were likely due to low early melt season surface 
albedos, measured by the Moderate Resolution Imaging Spectroradiometer (MODIS) and likely 
associated with the 31 March 2009 eruption of Mount Redoubt, southwestern Alaska. GRACE data 
compared well with in situ measurements at Wolverine Glacier (maritime Alaska), but poorly with those 
at Gulkana Glacier (interior Alaska). We conclude that, although GOA mass estimates from GRACE are 
robust over the entire domain, further constraints on subregional and seasonal estimates are necessary 
to improve fidelity to ground observations. 


INTRODUCTION 

Runoff from glaciers surrounding the Gulf of Alaska (GOA) 
impacts rates of global sea level (Gardner and others, 2013) 
and crustal uplift (Fu and Freymueller, 2012), and is the 
primary source of fresh water to terrestrial and aquatic 
ecosystems in the region (Neal and others, 2010). Unlike 
other geodetic mapping methods that only detect multi- 
annual changes (e.g. Arendt and others, 2002; Berthier and 
others, 2010), data from the NASA/German Aerospace 
Research Center (DLR) Gravity Recovery and Climate 
Experiment (GRACE) have the potential to resolve temporal 
variations in glacier mass, providing a unique tool for both 
annual and sub-annual mass-balance investigations. The 
primary challenges in working with GRACE data are their 
coarse spatial resolution and the contamination of the 
glacier gravity signal from other geophysical sources. To 
date, GRACE has been used primarily to quantify the 
cumulative mass balance of GOA glaciers in order to 
determine the total contribution of this region to rising sea 
level (Tamisiea and others, 2005; Chen and others, 2006a; 
Luthcke and others, 2008; Pritchard and others, 2010; Wu 
and others, 2010; Jacob and others, 2012; Sasgen and 
others, 2012a). Most studies have focused on the long-term 
trend in a GRACE time series because it is much less affected 
by errors in the modeling of terrestrial water storage (TWS). 
Most of the variability in TWS in the GOA region occurs at 


sub-annual timescales that match variations in glacier mass, 
making it more difficult to discriminate between glacio- 
logical and other sources of mass variation. 

Recent efforts have been made to compare GRACE GOA 
glacier mass-balance estimates with independent obser- 
vations. GOA mass-balance estimates derived from a high- 
resolution mass concentration (mascon) approach (Luthcke 
and others, 2008) have previously been validated against 
observations from aircraft altimetry (Arendt and others, 
2008) and weather station measurements of air temperature 
and precipitation (Arendt and others, 2009). Chen and 
others (2006b) compared a GRACE solution sampled at the 
location of two US Geological Survey (USGS) benchmark 
glaciers and found good agreement with field mass-balance 
observations. Other work (Hill and others, 201 1 ) has pointed 
toward potential errors in existing GRACE solutions for GOA 
glaciers, especially at sub-GOA basin scales where mass- 
balance magnitudes exceed what is expected from geodetic 
data (Berthier and others, 2010). As studies attempt to 
extract further spatial resolution from GRACE solutions (e.g. 
Jacob and others, 2012) there is an even greater need to 
constrain and validate higher-resolution mass-balance esti- 
mates for the GOA and other glacier regions. 

In this study we analyze the GOA glacier mascons from a 
new GRACE global mascon solution covering the period 
December 2003 to December 2010. This new solution is a 
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Fig- 1- GRACE solution domain for GOA glaciers. Boxes with gray outlines delineate 1 x 1 equal area degree mascons, and red outlines show 
groupings of these mascons into eight mountain ranges: 1. southwestern Alaska Range; 2. central Alaska Range; 3. eastern Alaska Range; 
4. Chugach Mountains; 5. Wrangell Mountains; 6. St Elias Mountains; 7. Juneau Icefield; 8. Stikine Icefield. Blue shading shows the location of 
Alaska glaciers, derived from a global glacier inventory (Arendt and others, 2012). Black dots show locations of ICESat elevation change 
measurements. Red dots show locations of USGS benchmark glaciers, and of Mount Redoubt volcano that erupted on 31 March 2009. 


subset of the high-resolution global mascon solution 
developed by the Space Geodesy Laboratory at the NASA 
Goddard Space Flight Center (Luthcke and others, 2013). 
We label this solution 'vl 2' following the version history of 
mascon solutions described in the literature (Luthcke and 
others, 2006a, 2008, 2013; Rowlands and others, 2010; 
Sabaka and others, 2010). Our goal is to compare temporal 
(annual to semi-annual) and spatial variations in the v12 
solution with independent data sources to assess the 
accuracy of the local mass-balance signal. We compare 
the v12 solution with a series of satellite altimetric, in situ 
mass-balance and climate datasets that serve both to 
validate GRACE and to provide insight into the geophysical 
drivers of seasonal and interannual GOA glacier mass 
balances during the GRACE observation period. 

DATA AND METHODOLOGY 
Solution domain 

We calculate gravity variations within a solution domain 
that encompasses most glaciers in Alaska, USA, as well as 
glaciers in Yukon and British Columbia, Canada, that are 
part of the icefields that straddle the US/Canada border (Fig. 
1 ). On its southeastern boundary, our domain ends just south 
of the southern end of the Stikine Icefield. We exclude those 
glaciers of the Aleutian Island Arc (1 793 km 2 ) and the Brooks 
Range of northern Alaska (531 km 2 ). The total glacier- 
covered area of our domain is 82 505 km 2 , calculated from 
a new glacier inventory (Randolph Glacier Inventory version 
2.0 or RGIv2.0) derived almost entirely from modern 
(~2005— 1 1 ) satellite imagery (Arendt and others, 2012). 
We solve for GRACE mascon parameters at every 1 x 1 
equal area mascon defined in Luthcke and others (2013). 
Our GOA region is a subset of the larger GOA domain 
defined in Luthcke and others (2013) that includes glaciers 
of the western Canadian Rockies. We exclude those glaciers 
in this analysis because they have traditionally not been 
included in earlier GOA mass-balance assessments. As a 
result, our mass-balance estimates are slightly different from 


those reported in Luthcke and others (2013). We combine 
mascons into eight mountain ranges, similar to those defined 
in Arendt and others (2002). 

GRACE solution 

The GRACE mission maps changes in Earth's gravity field by 
precisely measuring changes in distance between tandem 
satellites using a microwave K-band inter-satellite range and 
range-rate (KBRR) system with GPS data for positioning, 
star tracker data for orientation and accelerometer data 
to remove the non-conservative surface forces needed to 
isolate the gravity signal. GRACE measures changes in 
gravity from all components of the Earth system. In order 
to isolate the glacier mass-balance signal, it is necessary to 
forward-model the time-varying gravity from atmosphere, 
ocean, TWS and solid Earth variations. 

The v12 GRACE solution (Luthcke and others, 2013) 
calculates surface mass anomalies (mascons) directly from 
the GRACE KBRR observations taking into account the full 
noise covariance. We employ a unique method for accel- 
erometer calibration (Luthcke and others, 2006b) and for 
calculating mass changes as scale factors on the differential 
set of Stokes coefficients in a spherical harmonic expansion 
of the geopotential field (Rowlands and others, 2005). 
The mascons are estimated globally at 1 x 1 arcdeg 
(^25 000 km 2 ) spatial and 10 day temporal sampling, and 
we apply anisotropic temporal and spatial constraints 
between neighboring mascons representing similar surface 
types (land, ocean and glaciers). These neighbor constraints 
help to isolate signal and minimize signal leakage in and out 
of the land ice regions of interest. A detailed error and 
resolution analysis has shown that the basic mascon spatial 
resolution within a constraint region is equivalent to a 
Gaussian spatial smoother with 300 km radius (Luthcke and 
others, 2013). The signal within the GOA glacier region is 
isolated using these anisotropic spatial constraints, and the 
leakage errors estimated in Luthcke and others (2013) are 
included in our error estimates reported here. We estimate 
all errors at the 68% confidence (Icr) level. 
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The v12 mascon solution uses a similar GOA mascon 
solution domain and similar forward-modeling techniques 
to those in Luthcke and others (2008), including explicit 
accounting for post-Little Ice Age isostatic rebound resulting 
from the collapse of the Glacier Bay ice field (Larsen and 
others, 2005). A different ocean model is used, namely the 
Ocean Model for Circulation and Tides (Thomas, 2002), 
based on the Hamburg Ocean Primitive Equation Model 
(Drijfhout and others, 1996; Wolff and others, 1997). Over 
land we continue to use the 0.25° spatial, 3 hour temporal 
resolution GLDAS (Global Land Data Assimilation System) 
hydrology model (Rodell and others, 2004) with a 0.25° 
mask to remove glacierized regions from the TWS correc- 
tions, because these are not well modeled in GLDAS. We 
use a glacial isostatic adjustment (GIA) model based on the 
ICE-5G deglaciation history (Peltier, 2004), and an incom- 
pressible two-layer approximation to the VM2 viscosity 
profile (Paulson and others, 2007, computed and provided 
by Geruo A). We also apply geocenter corrections derived 
from the degree 1 Stokes coefficients determined from 
Swenson and others (2008). The v12 solution has been 
iterated three times to minimize residuals in modeled and 
observed KBRR data, resulting in improved recovery of non- 
modeled or mismodeled signal. 

We present 10 day mass anomaly estimates together with 
a 1 0day width Gaussian filter applied to the time series. We 
assume the Icr errors for each 10 day estimate are the 
difference between the filtered and unfiltered estimates, 
which we combine in a root-sum-square fashion to estimate 
the total error of the GOA mass changes due to noise in the 
GRACE signal. These errors are combined with leakage and 
forward-modeling errors described in Luthcke and others 
(2013) to obtain the total GOA error reported here. We 
calculate seasonal balances as the difference between yearly 
maximum and minimum values in the Gaussian-filtered 
mascon time series, and we identify balance years based on 
the end-of-summer minima. For mass-balance calculations 
over the entire period of record we report both an average of 
the annual balances determined by summation of the 
summer and winter balances, as well as mass trends 
determined from a least-squares simultaneous estimate of 
bias, trend, annual, semi-annual and 161 day cycle (alias 
period of S2 tide errors). 

We present summations of 1 x 1 mascon solutions over 
each of the eight mountain ranges, even though these 
mascon aggregates are not formally constrained in the way 
the total ice, land and ocean constraints are applied in the 
mascon approach (Luthcke and others, 2013). Our purpose 
is to ascertain the extent to which we can resolve mass 
variations at spatial scales smaller than the entire GOA 
region. Without any formal constraints, we are unable to 
quantify errors in the GRACE results we present for each of 
the eight regions. Instead, our comparison of the subregional 
GRACE estimates to independent climate and mass-balance 
observations is itself a form of error assessment, from which 
we can draw conclusions about the spatial resolution of the 
mascon approach. 

ICESat data 

We use Release 633 data from the GLA06 altimetry product 
of NASA's Geoscience Laser Altimeter System (GLAS) on 
board the Ice, Cloud and land Elevation Satellite (ICESat; 
Zwally and others, 2002) to calculate changes in surface 
elevation of GOA glaciers, and compare these with our 


GRACE-derived mass balances. We choose ICESat data for 
this purpose because they cover the entire GOA during 
2003-09, thereby fully overlapping with the GRACE obser- 
vation period. Other elevation datasets derived from air- 
borne altimetry (Johnson and others, 2013) or satellite 
stereoscopy (Berthier and others, 2010) are also available, 
but these span discontinuous time periods and do not cover 
the entire study region. The coverage of ICESat tracks over 
GOA glaciers is shown in Figure 1 . 

ICESat data were acquired over 1 7 repeated campaigns of 
a 33 day orbit sub-cycle between October 2003 and 
October 2009. ICESat provides surface elevations of 
~70m diameter laser altimetry footprints, each having a 
~170m along-track spacing. The elevation accuracy is 
highly dependent on surface slope but has been found to be 
within 1 m in typical glacier terrain (Moholdt and others, 
2010). In addition, there is a recently discovered bias in 
ICESat elevations resulting from an incorrect selection of the 
centroid, rather than the Gaussian peak, of the transmit 
pulse to calculate surface elevation (the GC offset; personal 
communication from A. Borsa, 2013). We simulate the 
effect of this time-variable offset using global mean offsets 
calculated for each observation campaign. 

We calculate elevation changes by examining elevation 
trends and residuals of surface planes fitted to 700 m long 
segments of near repeat-track data with 50% overlap, 
following methods in previous studies (Smith and others, 
2009; Moholdt and others, 2010, 2012; Gardner and others, 
2012, 2013). The method simultaneously solves for cross- 
track slope effects that are used to correct for repeat-track 
separation, as well as the true elevation changes used in our 
analysis. We use the same data-filtering parameters to 
remove data outliers as detailed in Moholdt and others 
(2010), except that we increase from 5 m to 10 m the 
maximum allowable elevation residual with respect to the 
trend of the plane. This was done to capture more of the 
steep topography and rapid elevation changes that are 
characteristic of many GOA glaciers. 

A total of 80 000 ICESat laser returns passed our filtering 
criteria, resulting in 6257 elevation change estimates that 
each cover a time-span of minimum 2 years within the 2003- 
09 period. We examined the sampling density within each of 
the eight mountain ranges, but found that the coverage of 
ICESat tracks was not sufficient for reliable mass-balance 
estimation in the smallest GOA regions. However, at the 
scale of the entire GOA region, we found that ICESat 
sampling density as a function of elevation matched well 
with the glacier area-elevation distribution (Fig. 2b), allow- 
ing us to proceed with regional volume change calculations. 
For each 200 m elevation bin we calculated differences 
between the ICESat sampling and area distribution density, 
and used that as a weighting factor on the ICESat elevation 
changes. This corrected elevation change was multiplied by 
the total glacier surface area within each bin to yield the total 
volume change in the GOA region. We converted the volume 
changes to mass changes assuming Sorge's law (Bader, 1 954), 
and an average density of 900 kg m -3 . We also calculated 
annual mass changes between ICESat campaigns in October/ 
November, which resulted in larger errors but allowed us to 
investigate links to our GRACE and climate observations. 

We calculated errors by dividing the standard deviation 
of elevation change by the square root of the number of 
uncorrelated observations, assuming a correlation length of 
5 km (Moholdt and others, 2012). We assume a fractional 
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Fig. 2. (a) Change in elevation determined from 2003-09 ICESat data, 
averaged within 200 m elevation bins. Error bars are the standard 
deviation of the mean change per elevation bin, divided by the 
square root of the number of uncorrelated observations, 
(b) Histograms of GOA glacier surface area (gray bars) and numbers 
of ICESat laser returns (orange line) used to estimate GOA glacier 
mass balance. The x-axis labels the midpoint of 200 m elevation bins 
over which the histograms and elevation changes are calculated. 


area uncertainty of ±10% (Arendt and others, 2006) and a 
density uncertainty of ±100kgrrr 3 (Moholdt and others, 
2012). The total mass-balance uncertainty was then calcu- 
lated from the quadrature sum of these error components. 
Sensitivity tests revealed that our error budget was domin- 
ated by the large scatter in ICESat elevation changes, with 
uncertainties in density and area having relatively small 
effects on the total error. Also, the uncertainty of the annual 
mass balance for 2009 was higher than for the other years 
due to the failure of the last laser instrument one-third of the 
way into the October 2009 observation campaign. 

Satellite snow-cover and albedo data 

We use the Moderate Resolution Imaging Spectoradiometer 
(MODIS)/Terra Snow Cover Monthly L3 Global MOD10CM 
V005 dataset (Hall and others, 2011) to investigate the 
spatial extent of snow cover on glacier-free terrain during the 
GRACE observation period. The MODIS snow-cover data 
are distributed at 0.05° spatial resolution of the Climate 
Modeling Grid. Daily maps are determined from a nor- 
malized difference snow index and averaged for the monthly 
product. The daily MODIS snow map accuracy has been 
reported at 93%, with larger errors occurring over complex 
terrain, and where snow cover is thin (Hall and Riggs, 2007). 
The largest source of errors likely results from clouds and 
mixed pixels (Painter and others, 2009), which tend to 
confuse the algorithm, generally resulting in an under- 
estimation of snow cover. 

We sample the snow-cover grids at locations within our 
GOA solution domain exclusive of the ice-covered regions. 
Our purpose is to assess interannual patterns in snow 
accumulation that may not be captured in ground station 
observations or reanalysis products, and to quantify poten- 
tial snow water equivalent errors in the GLDAS TWS model. 
We note that the 0.25° spatial, 3 hour temporal resolution 
version of the GLDAS model used here ingests the daily 
MODIS snow-cover product and adds an arbitrary snow 
water equivalent of 0.01 kgrrr 2 to locations at which 


MODIS gridcel Is show >40% snow cover (Rodell and 
Houser, 2004). Here we assess the sensitivity of TWS 
calculations to the 40% threshold value by counting the 
number of MODIS gridcel Is classified above and below 40% 
snow cover. 

We examine change in the spectral albedo of the glacier 
surfaces using the MODIS albedo product MCD43C4 for the 
years 2000-10. The MCD43C4 albedo is determined from 
the angular integration of a parametric model of the 
bidirectional reflectance distribution function (BRDF) that 
has been fit to 16 days of MODIS multi-angular reflectance 
observations in seven spectral bands. The product is 
available at 8 day intervals at a 0.05° horizontal resolution. 
We choose the coarse-resolution constant degree product 
over the higher-resolution 1000 m product (MCD43B3) for 
ease of processing, and because our GRACE analysis occurs 
on a low-resolution grid of mascons. The MODIS albedo 
product was found to agree with ground observations of 
snow albedo over the Greenland ice sheet within measure- 
ment error (Stroeve and others, 2005). However, particular 
attention must be given to the quality flags (QF) that 
accompany the dataset (Schaaf and others, 201 1). Here we 
only use albedo estimated from fully inverted BRDFs (QF = 
1). To ensure that we primarily capture changes in glacier 
albedo and not changes in terrestrial snow extent, we mask 
out all MODIS pixels with <95% glacier coverage as 
determined from the RGIv2.0. 

Climate data 

Here we build upon previous efforts to simulate seasonal 
variability of Alaska glaciers (Rasmussen and Conway, 2004; 
Arendt and others, 2009; Rasmussen and others, 2011), 
including an assessment of the relative influence of land 
surface physics on air temperature patterns, by comparing 
surface and 700 mbar air temperature records. For the upper 
air analysis we calculate the mean of two daily 700 mbar 
temperatures, one from the US National Centers for Environ- 
mental Prediction (NCEP) Reanalysis-2 and the second from 
the European Centre for Medium-Range Weather Forecasts 
(ECMWF) ERA-Interim reanalysis products. For each balance 
year in the GRACE record we calculate glacier area-weighted 
summer (June-August) departures from the 2004-10 mean 
for each of the eight mountain ranges. For the surface 
products we include both monthly temperature and precipi- 
tation fields derived from 1961-2009 station data using a 
total of 322 and 261 weather stations for temperature and 
precipitation respectively (note that 2010 data were not 
available). This surface product was calculated from absolute 
(for temperature) and proportional (for precipitation) anoma- 
lies of the station data from 1971-2000 climate normals 
(personal communication from D. Hill and S. Calos, 2011). 
These scattered anomalies were then interpolated onto a 
regular 2 km x 2 km grid covering all of Alaska. This 
approach assumes that the temporal derivative of the 
anomaly field is more spatially coherent than that of 
the normal field. The climate normals were taken from the 
Parameter-elevation Regressions on Independent Slopes 
Model (PRISM; Daly and others, 1994). PRISM uses a digital 
elevation model to correct for orographic effects on precipi- 
tation distribution, using regression equations that incorpor- 
ate elevation and aspect for each gridcell. We assume 
precipitation is a suitable proxy for snowfall, given the 
difficulties in determining a solid precipitation temperature 
threshold when using monthly data. Summer average 
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Fig. 3. Cumulative mass balance of GOA glaciers determined from 
the high-resolution GRACE v12 mascon solution (Luthcke and 
others, 2013). Ten-day estimates (blue dots), and the same data 
smoothed using a Gaussian filter with a 10 day averaging window 
(green curve); and trend (red line) recovered from simultaneous 
estimation of bias, trend, annual, semi-annual and 161 day cycle. 
GRACE errors for each 1 0 day estimate are assumed to be the differ- 
ence between the filtered (green line) and unfiltered (blue dot) data. 

(June-August) temperature and winter (September-May) total 
precipitation were calculated to compare with summer and 
winter mass balances determined from GRACE. We calcu- 
lated temperature and precipitation anomalies as the differ- 
ence of these average/totalized values from their mean values 
over the period of record. 

Field mass balance 

We compare GRACE measurements to field observations of 
glacier mass balance acquired by the USGS at Gulkana and 
Wolverine Glaciers (Van Beusekom and others, 2010). 
Biannual observations are carried out at each of three Index 
sites' that are used to represent different zones of each 
glacier. The glacier geometry and extent are updated every 
years. The glaciers are located in two distinctly different 
climate settings (Fig. 1). Gulkana Glacier had an area of 
15 km 2 in 2009 and is located on the south flank of the 
eastern Alaska Range in a continental climate. Wolverine 
Glacier was 17 km 2 in 2009 and is located in the Kenai 
Mountains in a wet, maritime climate. Measurements on 
both glaciers have been extrapolated to produce estimates of 
total GOA mass losses (Meier, 1984; Meier and Dyurgerov, 
2002) that are consistent with independent geodetic esti- 
mates (Cox and March, 2004; Harrison and others, 2009). 

RESULTS 

GRACE-derived mass balance 

Based on our least-squares estimate of the linear trend, the 
average annual mass balance (B a ) of GOA glaciers was 
—65 d= 1 1 Gta -1 during December 2003 to December 2010 
(Fig. 3). When we average only the 2004-10 B a values we 
obtain —71 ±11 Gta -1 . The difference results from the 
slightly different time period and the biases introduced in the 
linear fitting due to exceptionally negative summer balances 
at the start and end of the time series. We report both values 
because most GRACE studies calculate trends using the 
least-squares fitting method. 


Table 1. Annual mass balances (B a ) of GOA glaciers between 2004 
and 2010 balance years in units of mass (Gta -1 ), and in specific 
units (kgrrr 2 a -1 ) calculated by dividing by the total ice area in the 
mascon solution domain (82 505 km 2 ). Balance years begin during 
the fall of the previous calendar year. B a from GRACE is calculated 
as the difference between successive annual minima in Figure 3. B a 
from ICESat is derived from average elevation changes between fall 
campaigns in October/November each year 



GRACE 

ICESat 

Balance year 

B a 

B a 

B a B a 


Gt a -1 

kg m -2 a -1 

Gt a -1 kg rrr 2 a -1 

2004 

-1 25.0 ±12 

-1 .50 ± 0.1 5 

-11 4.0 ±31 -1 .4 0± 0.38 

2005 

—66.8 ±11 

-0.81 ±0.13 

-32.3 ±35 -0.39 ±0.43 

2006 

-54.5 ± 1 1 

-0.66 ±0.1 3 

-69.3 ±28 -0.85 ±0.34 

2007 

—66.5 ±11 

-0.81 ±0.13 

-11 5.0 ±33 -1.40 ±0.40 

2008 

1 5.0 ± 1 0 

0.1 8 ± 0.12 

16.5 ±39 0.20 ±0.47 

2009 

-1 59.0 ±28 

-1 .9 0± 0.34 

-11 9.0 ±99 -1.40 ±1.20 

2010 

-43.1 ±28 

-0.52 ±0.34 


2004-09* 

—76.0 ±11 

-0.92 ±0.1 3 

-72.0 ±22 -0.87 ± 0.26 

2004-09^ 

—61 .0 ± 1 1 

-0.74 ±0.1 3 

-65.0 ±12 -0.79 ± 0.1 5 

2004-10* 

—71 .0 ± 1 1 

-0.86 ±0.1 3 


2004-1 Qt 

-65.0 ± 1 1 

-0.78 ±0.1 3 



* Average of seasonal values. ^Trend using all values. 


The range in B a in Alaska between 2004 and 2010 was 
174Gt, approximately double the mass loss rate estimated 
from a linear fit to the entire time series (Table 1 ). B a became 
progressively less negative from 2004 to 2008, with near- 
balance conditions occurring in 2008. This was followed in 
2009 by the most negative B a in the 2004-10 period. We 
examined departures of summer (B s ) and winter (B w ) 
balances from the 2004-10 mean, calculated as the ratio of 
each year's seasonal balance to the mean seasonal balance 
over the period of the GRACE record. Our sign convention is 
such that a negative departure in each season indicates a 
more negative mass balance. B s departures (-83 to 91 Gta -1 ) 
were larger than B w departures (-25 to 22 Gta -1 ), showing 
that variations in B s accounted for most of the variability in B a 
(Fig. 4). The slightly positive 2008 and large negative 2009 B a 
values were caused by strong positive and negative 
departures of B s in the respective years. 

We observed large spatial variability in the patterns of 
glacier mass loss between different balance years (Fig. 5). 
Notable mass gains occurred for glaciers in the central and 
eastern Alaska ranges (mountain ranges 2 and 3) during the 
2008 balance year, the Stikine Icefield (mountain range 8) 
during 2007, and glaciers of the Alaska and Kenai 
peninsulas (mountain ranges 1 and 4) during 2008 (Fig. 1). 
Large mass losses occurred through all regions during the 
2004 and 2009 balance years due to exceptionally negative 
B s . During the entire period, the St Elias, Glacier Bay and 
Juneau icefields regions lost the greatest amount of mass 
relative to other regions. The amplitudes in glacier mass 
increased with decreased latitude, and were largest in 
southeast Alaska. 

Comparison to ICESat-derived mass balance 

ICESat elevation changes averaged over the 2003-09 
observation period show a strong elevation dependence, 
with 2-4 m a -1 of thinning at low elevations tapering to 
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Fig. 4. Departures (difference from 2004-1 0 mean) of GOA summer 
(orange) and winter (gray) mass balances. Seasonal balances are 
calculated as differences between successive minima and maxima 
in the GOA mascon time series (Fig. 3). Negative departures for 
both summer and winter balances indicate a contribution toward 
more negative annual balances, as compared to the 7 year mean. 


near-zero changes at high elevations (Fig. 2a). We multiplied 
these elevation changes by the area distribution of GOA 
glaciers to obtain an average B a of -65±12Gta _1 or 
— 0.79 ± 0.1 5 kg m -2 a -1 (Table 1). This is 4Gta _1 more 
negative than GRACE data subsampled to the same 
period as ICESat. GRACE and ICESat also agree to within 
error bars for individual balance years (B a ) except for 2007 


Fig. 6. Comparison of area-averaged annual mass balances in the 
GOA region derived from GRACE (orange) and ICESat (gray). 
Balance years begin in the fall of the previous calendar year, at the 
time of minimum mass (GRACE), or in October/November (ICESat). 


(Table 1; Fig. 6). Our simulations of the ICESat GC offset 
indicated that the total impact on the GOA mass balance 
was <1 Gta -1 . Because this was well within our calculated 
ICESat error estimates, we do not include this correction in 
our analysis. 

Comparison to field observations 

We compare field observations of B w , B s and B a at Gulkana 
and Wolverine Glaciers to their closest 1 x 1 mascon 



Year 


Winter 

Annual 

Summer 


r I | 


50 Gt 


Fig. 5. Spatial distribution of summer, winter and annual mass balance (orange, gray bars and red lines) in Alaska from the GRACE v12 
mascon solution (Luthcke and others, 2013), for each balance year spanning 2004-10. Mountain ranges (numbered 1-8) are labeled in 
Figure 1. Summer, winter and annual balances are calculated as in Table 1. Blue shading shows glacier locations. 
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Fig. 7. Comparison of area-averaged summer (orange squares), winter (gray circles) and annual (red diamonds) mass balances from GRACE and 
those derived from USGS field observations at (a) Gulkana Glacier and (b) Wolverine Glacier between 2003 and 2009. GRACE data are 
converted to specific values by dividing the mass change by the total ice area in mascons containing each glacier. 


(Fig. 7). We convert the GRACE mass changes to specific 
values by dividing by the total ice area in each mascon. Our 
scatter plots illustrate both the strength of the correlation 
between the two datasets (i.e. similarity of patterns from one 
year to the next) as well as the level of agreement in 
magnitude between GRACE and ground observations (meas- 
ured by the departure of points from the diagonal one-to-one 
line). At Gulkana Glacier, B a from GRACE and ground 
observations were poorly correlated (r 2 = 0.05, p = 0.1), 
and the absolute values of GRACE B s and B w values were 
larger than ground observations. For Wolverine Glacier, B a 
had a relatively high correlation (r 2 = 0.75, p = 0.79) 
between ground and GRACE observations, and the magni- 
tudes of field-derived B s and B w were very close to those 
observed by GRACE. We note that the higher correlation 
between GRACE- and field-derived B a values for Wolverine 
Glacier may have occurred due to compensating errors in 
the associated B s and B w values, as indicated by their lower 
correlation to the GRACE data. 

Comparisons to climate observations 

Both the 700mbar reanalysis air temperatures and the 
interpolated ground surface temperatures captured a large 
amount of the variability in GRACE-derived B s data 
(r 2 as 0.72, p = 0.02 and r 2 — 0.59, p = 0.07 respectively; 
Fig. 8; note the inverted secondary y-axis). In nearly all cases, 
surface and 700mbar anomalies were positive during years 
of negative B s . A notable exception was 2009, a year of 
strongly negative B s during which 700mbar temperatures 
were slightly above normal, but surface temperatures were 
slightly below normal. This was in contrast to 2004, when a 
similarly large negative B s value was matched by strongly 
positive surface and 700mbar anomalies. We note that 
during 2004 the magnitude of the 700mbar temperature 
anomaly was approximately four times larger than in 2009, 
even though B s anomalies were similar between the two 
years. During 2004, 700mbar air temperature departures 
from the 2004-10 mean were distributed across the GOA 
region, while in 2009 the highest departures occurred over 
the St Elias, Juneau and Stikine icefield regions (Fig. 9). 2008 
had the largest positive summer balance departure that was 
matched by the most negative surface and 700 mbar air tem- 
perature departures. The 2008 temperature departures were 
distributed relatively uniformly across the region (Fig. 9). 


There was no significant correlation between GRACE- 
derived B w and winter (September-May) precipitation 
determined from PRISM (r 2 =0.16, p = 0.43) and from 
the reanalysis products (r 2 = 0.07, p = 0.57; Fig. 10). The 
winter balance years 2006 and 2008 had above-normal 
winter balances together with above-normal surface snow 
accumulation values, while the opposite occurred for the 
2007 and 2009 balance years. 

Spectral albedo values averaged over May-July for GOA 
glacier surfaces showed a 50% increase in shortwave 
absorption during 2009 relative to the average of the 
remaining years in 2004-10 (Fig. 1 1). The greatest decrease 
in albedo occurred in the visible and near-infrared spectrum 
(0.40-0.9 pm). 

Snow cover on glacier-free terrain (an annual average of 
all snow-cover values in the GOA region exclusive of 
glaciers) had maximum values of 90-96% and minimum 
values of 7-13% within the 2004-10 time period (Fig. 12). 
The summer of 2008 had the largest minimum snow cover 
(13%), corresponding with the least negative B s discussed 



Fig. 8. Departures of GRACE GOA summer glacier mass balances 
from the 2004-10 mean (orange circles, left axis; note that negative 
summer balance departures imply greater summer mass loss). GOA 
summer (June-August) surface (black squares) and 700 mbar lower- 
troposphere (black circles) air temperature departures from the 
2004-09 mean (right axis; note inverted scale). 
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Fig. 9. Departures in summer (June-August) 700 mbar air temperature (mean of NCEP R-2 and ERA-Interim reanalysis) from 2004-1 0 mean at 
each of the eight GRACE mascon mountain ranges. 


above. For most years the average snow cover was <40%, 
and therefore below the threshold for inclusion in the GLDAS 
model, during May-July. To assess potential errors in 
eliminating cells with <40% snow cover, we multiplied the 
number of cells with <40% snow cover by 0.01 kgrrr 2 
(Rodell and EHouser, 2004), yielding a value of 2.5 Gt. We 
lack further information to more accurately assess the water 
equivalent depth of snow remaining, but we provide this 


estimate to illustrate the potential order-of-magnitude impact 
of late-season snowpacks on our mass trend retrievals. 

DISCUSSION 

Our mass-balance estimates averaged over the 2004-09 
balance years from GRACE and ICESat compare well, 
providing mutual validation of the magnitude of GOA glacier 
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Fig. 10. Departures of GRACE GOA winter mass balances from the 
2004-10 mean (gray circles, left axis). GOA winter (September- 
May) total winter snowfall departure from the 2004-09 mean (right 
axis). Data are derived from station observations (black squares) 
and from reanalysis model output (red squares). 


contributions to rising sea level during this period. Our 
estimates are more negative than two other GRACE estimates 
for GOA glaciers over the same time-span (—42 =t 6Gta _1 
(Jacob and others, 2012) and -54±13Gta _1 (Sasgen and 
others, 2012b)). The relatively strong agreement between 
GRACE- and ICESat-derived B a for individual years provides 
good support for our assessment that interannual variations in 
GOA GRACE signals are due to fluctuations in glacier mass. 
The elevation dependence of ICESat elevation changes also 
matches recent airborne altimetry observations (e.g. Johnson 
and others, 201 3), lending further support to the accuracy of 
our ICESat analysis for GOA glaciers. 

Both surface and ZOOmbar air temperatures explain 
much of the variability in B s as observed by GRACE. The 
slightly stronger correlation for 700mbar suggests lower- 
troposphere data may be more suitable than surface data for 
analysis of snow and ice variations, supporting previous 
studies (Rasmussen and Conway, 2004). This is likely 
because 700 mbar temperatures are less subject to localized 
surface effects and are more indicative of regional tempera- 
ture variations. The relatively weak correlation of £> w with 
snow accumulation observations likely results from diffi- 
culties in quantifying snow cover from precipitation data, 
and the paucity of ground observations to constrain the 
reanalysis and surface precipitation products. 

We attribute the strongly negative B s in 2009 to the 
31 March 2009 eruption of Mount Redoubt that spread ash 
over a large area, causing a significant reduction in surface 
albedo as measured by MODIS. Ash dispersal maps from the 
eruption show the plume progressing primarily north and 
east of the eruption center, with ash fall covering an 
estimated 80 000 km 2 in the vicinity of Mount Redoubt 
(Wallace and others, 2013). Trace ash (particle size 
<0.8 mm) was observed as far away as Fairbanks, Alaska, 
located 550 km north-northeast of the volcano (Wallace and 
others, 2013). These observations are consistent with our 
2009 field operations, during which several of us observed 
ash as far away as the central Alaska and St Elias Ranges 
(^300 and 650 km from Mount Redoubt, respectively). 
Snow albedo is very sensitive to the presence of highly 
absorbing atmospheric aerosols such as those expelled 



Fig. 11. Average May-July MCD43C3 MODIS-derived spectral 
albedo for all gridcel Is within a 250km radius of Mount Redoubt 
and having a glacier coverage of 95% or greater for the years 
2000-10. Areas listed next to sample years indicate the total 
gridcell area sampled each year. 

during a volcanic event (Conway and others, 1996). We 
suggest the Redoubt eruption caused an enhancement of 
surface melt through greater absorption of solar radiation, 
providing a mechanism for increased GOA glacier mass loss 
even in the absence of higher air temperatures. 

The sequence and timing of air temperature and snow 
accumulation anomalies appear to be important in control- 
ling GOA mass-balance patterns. The relatively large winter 
accumulation during the 2008 balance year was followed by 
a cool summer that allowed snow to remain unmelted on 
glacier-free surfaces. The fact that most of the GOA mascons 
had <40% snow-cover values during the summer months 
means that they were not accounted for in the GLDAS 
model (Rodell and Houser, 2004). In the event that summer 



Fig. 12. Average percent snow cover on glacier-free terrain in the 
mascon solution domain determined from the MODIS Monthly L3 
Global snow-cover dataset. Red line is the monthly time series, and 
gray bars mark the magnitude of the summer minimum value. 
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snowpacks registering <40% cover in the MODIS product 
are of sufficient extent or thickness to be sensed by the 
GRACE satellites, this is a bias in our calculations that 
should be addressed in future work. 

B s and B w values predicted by GRACE are much larger 
than is physically realistic for Gulkana Glacier, and far 
exceed any in situ seasonal balances measured during its 
entire record dating back to 1966. These erroneous GRACE 
amplitudes likely result from some combination of signal 
leakage across the solution domain (Sabaka and others, 
2010) and signal contamination due to TWS variations on 
glacier-free surfaces. Concerning signal leakage, we note 
that although our GRACE solution does have spatial and 
temporal constraints applied to the entire GOA solution 
domain, no constraints have been applied to either the 
individual mascons or mountain ranges. These constraints 
help to minimize leakage of signal between different broad- 
scale geophysical systems, such as land, glacier and ocean 
surface types. Without constraints, individual mascons can 
contain signal that is to some extent smeared across a larger 
region of the Earth's surface. For example, Luthcke and 
others (2013) found that signal from a single mascon can 
affect other mascons at spatial scales up to 600 km. These 
problems reflect fundamental limitations of the GRACE 
satellites that are made worse when applying GRACE to a 
complex region such as the GOA that has a highly non- 
uniform distribution of glaciers near a highly dynamic 
ocean boundary. The fact that our comparisons with 
Gulkana Glacier, located in a region with sparse glacieriza- 
tion, were much worse than with Wolverine Glacier, 
located in a heavily glacierized region, suggests that signals 
from the large ice masses in the St Elias and Chugach 
Mountains are leaking into mascons in the Alaska Range 
and other regions. 

TWS modeling errors could also account for discrep- 
ancies between ground observations and GRACE. Luthcke 
and others (2008) showed that forward modeling of TWS 
using the same dataset employed for our v12 solutions (the 
GLDAS/NOAH product) did not match GRACE observations 
over much of the glacier-free surfaces of interior Alaska 
during 2003-08. TWS errors result from a paucity of snow- 
depth, groundwater and other observations needed to 
calibrate the GLDAS product. If TWS is mismodeled over 
glacier-free surfaces located within glacier mascons, this will 
result in erroneous attribution of terrestrial snow cover, 
groundwater and other TWS signals to glacier surfaces. Such 
errors have the potential to be greatest for glacier mascons 
with low concentrations of glacier ice, wherein a greater 
proportion of TWS corrections are being applied. We note 
that for those years when all snow cover melts on glacier- 
free terrain, these errors only affect the sub-annual com- 
ponents of mass change, and not the annual trend. 


CONCLUSIONS 

GOA glaciers are controlled by numerous environmental 
forcings that complicate interpretation and modeling of their 
mass changes, especially over short time periods. During 
2003-08, the mass balance of GOA glaciers was strongly 
correlated to mean summer air temperatures, emphasizing 
the importance of acquiring air temperature data from high 
mountain regions in order to improve model predictions of 
Alaska glacier surface mass balances. At the same time, our 
2009 observations highlight potential limitations in models 


that only include temperature as a driver of glacier ablation. 
Due to the regional reduction in albedo from the 31 March 
2009 eruption of Mount Redoubt, temperature alone was an 
insufficient proxy for glacier surface mass balance, and other 
factors (e.g. changes to the optical properties of the snow 
and ice surface) may have been the dominant controls. 

Studies aimed at matching GRACE observations to field 
observations require consideration of leakage errors and 
fundamental resolution constraints inherent in satellite 
gravimetry datasets. We have shown that individual mascon 
time series combined additively for the purpose of hydro- 
logical analysis within mountain ranges generally do not 
match the magnitude of ground observations, even though 
they are highly temporally correlated. We attribute this to 
the smearing of signal within the GOA domain and the 
mismodeling of TWS variations, both of which are enhanced 
by the complex spatial distribution of glaciers within the 
region and the lack of sub-GOA solution constraints. 

Snow on ground likely plays an important role in 
controlling the mass balance of individual glacier catch- 
ments, and biases to GRACE-derived mass balances may 
occur when the summer snowpack on glacier-free terrain 
does not completely melt in a given year. Further improve- 
ments to TWS models and to the acquisition of accurate 
snow water equivalence data are important to improve 
recovery and attribution of GRACE mass-balance signals in 
this region. 

We have provided the first independent validation of 
GRACE glacier mass-balance estimates spanning the entire 
GOA region, using data from the ICESat mission. The 
agreement of these two estimates strongly supports our 
assertion that inaccuracies in sub-annual mass balances, due 
to problems outlined above, are sufficiently random over the 
length of the GRACE record that they do not bias our multi- 
annual estimates. This suggests that our existing GRACE 
processing constraints used to separate mass change signals 
occurring within the entire GOA region from surrounding 
land and ocean surfaces are robust. The relatively strong 
correlation between B a values derived from GRACE and 
ICESat indicates that future satellite altimetric missions, 
which will have an even higher sampling density than 
ICESat, will be a valuable tool for annual GOA mass-balance 
assessments. 
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